@testset "Instances" begin
  @testset "Perfect bipartite matching" begin
    @testset "Constructor with $i nodes on each side" for i in [2, 5, 10] 
      n = i
      epsilon = (i - 2) / 16
      reward = Distribution[Bernoulli(.5 + ((i == j) ? epsilon : 0.)) for i in 1:n, j in 1:n]
      instance = PerfectBipartiteMatching(reward, PerfectBipartiteMatchingNoSolver())

      @test instance.n_arms == n ^ 2
      @test instance.reward == reward
      @test optimal_reward(instance) == instance.optimal_reward
      
      # Error: non bipartite. 
      reward = Distribution[Bernoulli(.5 + ((i == j) ? epsilon : 0.)) for i in 1:n, j in 1:(n + 2)]
      @test_throws ErrorException PerfectBipartiteMatching(reward, PerfectBipartiteMatchingNoSolver())
      
      reward = Distribution[Bernoulli(.5 + ((i == j) ? epsilon : 0.)) for i in 1:(n + 2), j in 1:n]
      @test_throws ErrorException PerfectBipartiteMatching(reward, PerfectBipartiteMatchingNoSolver())
    end

    @testset "State with $i nodes on each side" for i in [2, 5, 10] 
      n = i
      epsilon = (i - 2) / 16
      reward = Distribution[Bernoulli(.5 + ((i == j) ? epsilon : .0)) for i in 1:n, j in 1:n]
      instance = PerfectBipartiteMatching(reward, PerfectBipartiteMatchingNoSolver())

      state = initial_state(instance)

      @test state.round == 0
      @test state.regret == 0.0
      @test state.reward == 0.0
      @test length(state.arm_counts) == n * n
      @test length(state.arm_reward) == n * n
      @test length(state.arm_average_reward) == n * n

      for i in 1:n
        for j in 1:n
          @test state.arm_counts[(i, j)] == 0
          @test state.arm_reward[(i, j)] == 0.0
          @test state.arm_average_reward[(i, j)] == 0.0
        end
      end
    end

    @testset "Trace with $i nodes on each side" for i in [2, 5, 10] 
      n = i
      epsilon = (i - 2) / 16
      reward = Distribution[Bernoulli(.5 + ((i == j) ? epsilon : .0)) for i in 1:n, j in 1:n]
      instance = PerfectBipartiteMatching(reward, PerfectBipartiteMatchingNoSolver())

      trace = initial_trace(instance)

      @test length(trace.states) == 0
      @test length(trace.arms) == 0
      @test length(trace.reward) == 0

      @test eltype(trace.states) == State{Tuple{Int, Int}}
      @test eltype(trace.arms) == Vector{Tuple{Int, Int}}
      @test eltype(trace.reward) == Vector{Float64}
    end

    @testset "Pull with $i nodes on each side" for i in [2, 5, 10]
      n = i
      reward = Distribution[Bernoulli(((i == j) ? 0.0 : 1.0)) for i in 1:n, j in 1:n]
      instance = PerfectBipartiteMatching(reward, PerfectBipartiteMatchingNoSolver())

      srand(1)
      @test pull(instance, [(1, 2), (i, i)]) == [1.0, 0.0]
    end

    @testset "Check feasibility with 3 nodes on each side" begin
      n = 3
      reward = Distribution[Bernoulli(((i == j) ? 0.0 : 1.0)) for i in 1:n, j in 1:n]
      instance = PerfectBipartiteMatching(reward, PerfectBipartiteMatchingNoSolver())

      @test is_feasible(instance, [(1, 2), (2, 3), (3, 1)])
      @test ! is_feasible(instance, [(2, 3), (3, 2), (1, 2)])
      @test ! is_feasible(instance, [(2, 3), (3, 2), (2, 1), (1, 2)])
    end

    @testset "LP solver" begin
      @testset "Constructor" for i in [2, 5, 10]
        n = i
        reward = Distribution[Bernoulli(((i == j) ? 0.0 : 1.0)) for i in 1:n, j in 1:n]
        instance = PerfectBipartiteMatching(reward, PerfectBipartiteMatchingLPSolver(CbcSolver(logLevel=0)))

        @test instance.solver != nothing
        @test instance.solver.m != nothing
        @test size(instance.solver.x, 1) == n
        @test size(instance.solver.x, 2) == n
      end

      @testset "Optimal reward with $i nodes on each side" for i in [2, 5, 10] 
        n = i
        epsilon = (i - 2) / 16
        reward = Distribution[Bernoulli(.5 + ((i == j) ? epsilon : 0.)) for i in 1:n, j in 1:n]
        instance = PerfectBipartiteMatching(reward, PerfectBipartiteMatchingLPSolver(CbcSolver(logLevel=0)))
  
        @test instance.optimal_reward == i * (.5 + epsilon) # Each time, choose the edge with mean .5 + epsilon. 
        @test optimal_reward(instance) == instance.optimal_reward
      end

      @testset "Solve with $i nodes on each side" for i in [2, 5, 10]
        n = i
        reward = Distribution[Bernoulli(((i == j) ? 0.0 : 1.0)) for i in 1:n, j in 1:n]
        instance = PerfectBipartiteMatching(reward, PerfectBipartiteMatchingLPSolver(CbcSolver(logLevel=0)))

        srand(i)
        drawn = Dict((i, j) => rand() for i in 1:n, j in 1:n)
        solution = solve_linear(instance, drawn)
        @test is_feasible(instance, solution)
      end
    end

    @testset "Munkres solver" begin
      @testset "Constructor" for i in [2, 5, 10]
        n = i
        reward = Distribution[Bernoulli(((i == j) ? 0.0 : 1.0)) for i in 1:n, j in 1:n]
        instance = PerfectBipartiteMatching(reward, PerfectBipartiteMatchingMunkresSolver())

        @test instance.solver != nothing
      end
      
      @testset "Optimal reward with $i nodes on each side" for i in [2, 5, 10] 
        n = i
        epsilon = (i - 2) / 16
        reward = Distribution[Bernoulli(.5 + ((i == j) ? epsilon : 0.)) for i in 1:n, j in 1:n]
        instance = PerfectBipartiteMatching(reward, PerfectBipartiteMatchingMunkresSolver())
  
        @test instance.optimal_reward == i * (.5 + epsilon) # Each time, choose the edge with mean .5 + epsilon. 
        @test optimal_reward(instance) == instance.optimal_reward
      end

      @testset "Solve with $i nodes on each side" for i in [2, 5, 10]
        n = i
        reward = Distribution[Bernoulli(((i == j) ? 0.0 : 1.0)) for i in 1:n, j in 1:n]
        instance = PerfectBipartiteMatching(reward, PerfectBipartiteMatchingMunkresSolver())

        srand(i)
        drawn = Dict((i, j) => rand() for i in 1:n, j in 1:n)
        solution = solve_linear(instance, drawn)
        @test is_feasible(instance, solution)
      end
    end

    @testset "Hungarian solver" begin
      @testset "Constructor" for i in [2, 5, 10]
        n = i
        reward = Distribution[Bernoulli(((i == j) ? 0.0 : 1.0)) for i in 1:n, j in 1:n]
        instance = PerfectBipartiteMatching(reward, PerfectBipartiteMatchingMunkresSolver())

        @test instance.solver != nothing
      end
      
      @testset "Optimal reward with $i nodes on each side" for i in [2, 5, 10] 
        n = i
        epsilon = (i - 2) / 16
        reward = Distribution[Bernoulli(.5 + ((i == j) ? epsilon : 0.)) for i in 1:n, j in 1:n]
        instance = PerfectBipartiteMatching(reward, PerfectBipartiteMatchingMunkresSolver())
  
        @test instance.optimal_reward == i * (.5 + epsilon) # Each time, choose the edge with mean .5 + epsilon. 
        @test optimal_reward(instance) == instance.optimal_reward
      end

      @testset "Solve with $i nodes on each side" for i in [2, 5, 10]
        n = i
        reward = Distribution[Bernoulli(((i == j) ? 0.0 : 1.0)) for i in 1:n, j in 1:n]
        instance = PerfectBipartiteMatching(reward, PerfectBipartiteMatchingMunkresSolver())

        srand(i)
        drawn = Dict((i, j) => rand() for i in 1:n, j in 1:n)
        solution = solve_linear(instance, drawn)
        @test is_feasible(instance, solution)
      end
    end
  end

  @testset "Shortest paths in a complete graph" begin
    @testset "Constructor with $i nodes" for i in [2, 5, 10] 
      # Edge weights: high gain likelihood for the path 1->2->3...n, low for the other edges. 
      n = i
      epsilon = 1 / n
      reward = Distribution[Bernoulli((i + 1 == j) ? (1 - epsilon) : (epsilon)) for i in 1:n, j in 1:n]
      instance = CompleteGraphShortestPath(reward, 1, n, CompleteGraphShortestPathNoSolver())

      @test instance.n_arms == n * (n - 1) / 2
      @test instance.reward == reward
      @test optimal_reward(instance) == instance.optimal_reward
      
      # Error: source makes no sense. 
      @test_throws ErrorException CompleteGraphShortestPath(reward, - 1, n, CompleteGraphShortestPathNoSolver())
      @test_throws ErrorException CompleteGraphShortestPath(reward, n + 1, n, CompleteGraphShortestPathNoSolver())
      
      # Error: destination makes no sense. 
      @test_throws ErrorException CompleteGraphShortestPath(reward, 1, - 1, CompleteGraphShortestPathNoSolver())
      @test_throws ErrorException CompleteGraphShortestPath(reward, 1, n + 1, CompleteGraphShortestPathNoSolver())
    end

    @testset "State with $i nodes" for i in [2, 5, 10] 
      n = i
      epsilon = 1 / n
      reward = Distribution[Bernoulli((i + 1 == j) ? (1 - epsilon) : (epsilon)) for i in 1:n, j in 1:n]
      instance = CompleteGraphShortestPath(reward, 1, n, CompleteGraphShortestPathNoSolver())

      state = initial_state(instance)

      @test state.round == 0
      @test state.regret == 0.0
      @test state.reward == 0.0
      @test length(state.arm_counts) == n * n
      @test length(state.arm_reward) == n * n
      @test length(state.arm_average_reward) == n * n

      for i in 1:n
        for j in 1:n
          @test state.arm_counts[(i, j)] == 0
          @test state.arm_reward[(i, j)] == 0.0
          @test state.arm_average_reward[(i, j)] == 0.0
        end
      end
    end

    @testset "Trace with $i nodes" for i in [2, 5, 10] 
      n = i
      epsilon = 1 / n
      reward = Distribution[Bernoulli((i + 1 == j) ? (1 - epsilon) : (epsilon)) for i in 1:n, j in 1:n]
      instance = CompleteGraphShortestPath(reward, 1, n, CompleteGraphShortestPathNoSolver())

      trace = initial_trace(instance)

      @test length(trace.states) == 0
      @test length(trace.arms) == 0
      @test length(trace.reward) == 0

      @test eltype(trace.states) == State{Tuple{Int, Int}}
      @test eltype(trace.arms) == Vector{Tuple{Int, Int}}
      @test eltype(trace.reward) == Vector{Float64}
    end

    @testset "Pull with $i nodes" for i in [2, 5, 10]
      n = i
      epsilon = 1 / n
      reward = Distribution[Bernoulli((i + 1 == j) ? (1 - epsilon) : (epsilon)) for i in 1:n, j in 1:n]
      instance = CompleteGraphShortestPath(reward, 1, n, CompleteGraphShortestPathNoSolver())

      srand(1)
      if i == 2
        @test pull(instance, [(1, 2), (i, i)]) == [1.0, 1.0]
      else
        @test pull(instance, [(1, 2), (i, i)]) == [1.0, 0.0]
      end
    end

    @testset "Check feasibility with 3 nodes" begin
      n = 3
      epsilon = 1 / n
      reward = Distribution[Bernoulli((i + 1 == j) ? (1 - epsilon) : (epsilon)) for i in 1:n, j in 1:n]
      instance = CompleteGraphShortestPath(reward, 1, n, CompleteGraphShortestPathNoSolver())

      @test is_feasible(instance, [(1, 3)])
      @test is_feasible(instance, [(1, 2), (2, 3)])
      @test ! is_feasible(instance, Tuple{Int, Int}[])
      @test ! is_feasible(instance, [(1, 2), (2, 3), (3, 1)])
      @test ! is_feasible(instance, [(2, 3), (3, 2), (1, 2)])
      @test ! is_feasible(instance, [(2, 3), (3, 2), (2, 1), (1, 2)])
    end

    @testset "LightGraphs.jl Dijkstra solver" begin
      @testset "Constructor with $i nodes" for i in [2, 5, 10]
        n = i
        epsilon = 1 / n
        reward = Distribution[Bernoulli((i + 1 == j) ? (1 - epsilon) : (epsilon)) for i in 1:n, j in 1:n]
        instance = CompleteGraphShortestPath(reward, 1, n, CompleteGraphShortestPathLightGraphsDijkstraSolver())

        @test instance.solver != nothing
        @test instance.solver.graph != nothing
        @test size(instance.solver.weights_matrix, 1) == n
        @test size(instance.solver.weights_matrix, 2) == n
        @test instance.solver.source == 1
        @test instance.solver.destination == n
      end

      @testset "Optimal reward with $i nodes" for i in [2, 5, 10] 
        n = i
        epsilon = 1 / n
        reward = Distribution[Bernoulli((i + 1 == j) ? (1 - epsilon) : (epsilon)) for i in 1:n, j in 1:n]
        instance = CompleteGraphShortestPath(reward, 1, n, CompleteGraphShortestPathLightGraphsDijkstraSolver())
  
        @test instance.optimal_reward ≈ (n - 1) * (1 - epsilon) atol=1.e-9 # Choose the path along the (1 - epsilon) edges. 
        @test optimal_reward(instance) == instance.optimal_reward
      end

      @testset "Solve with $i nodes" for i in [2, 5, 10]
        n = i
        epsilon = 1 / n
        reward = Distribution[Bernoulli((i + 1 == j) ? (1 - epsilon) : (epsilon)) for i in 1:n, j in 1:n]
        instance = CompleteGraphShortestPath(reward, 1, n, CompleteGraphShortestPathLightGraphsDijkstraSolver())

        srand(i)
        drawn = Dict((i, j) => rand() for i in 1:n, j in 1:n)
        solution = solve_linear(instance, drawn)
        @test is_feasible(instance, solution)
      end
    end
  end
end